library(ggplot2)
library(mnormt)

f<-function(z,y,a,Sig)
{
	x<-cbind(z,y)
	res=0
	res<-dmnorm(x,m=c(0,0),varcov=diag(2),log=TRUE)
	L<-chol(A)
	res[L[1,]%*%t(x)<a[1] | L[2,]%*%t(x)<a[2]]<- -Inf
	#res[x[,1]<a[1] | x[,2]<a[2]]<- -Inf
	return(res)
}
Sig<-A
dyn.load("test/Orthant.so")
b<-a
L<-chol(A)
res<-.C("EPtrunc",as.integer(2),as.double(t(L)),as.double(b),b=as.double(b),A=as.double(A))

x<-rmvnorm(5000,mean=as.matrix(res$b),sigma=matrix(res$A,ncol=2))
X <- seq(from = min(x[,1]), to = max(x[,1]), length.out = 100)
Y <- seq(from = min(x[,2]), to =max(x[,2]), length.out = 100)
Z <- outer(X, Y,f,a=a,Sig=Sig)
df<- expand.grid(X, Y)
datap<-as.data.frame(x)

df$logpost <- as.vector(Z)
names(df) <- c("mu1", "mu2", "logposterior")
g <- ggplot(df, aes(x = mu1, y = mu2)) + geom_tile(aes(fill = exp(logposterior)))
g <- g + stat_contour(aes(z = exp(logposterior)))
g <- g + xlab(expression(x[1])) + ylab(expression(x[2]))
g<-g+geom_point(data=datap,aes(datap$V1,datap$V2),colour="red",size=1)
